% Figure 
load '../data/DataforMatlab/UK_Budget.mat'
load '../data/DataforMatlab/pvec_tmp.mat'
subset = [1:9 11:16];
Bx_vec = sum(B_vec(subset,:,:),1);
Bxi_vec = B_vec(subset,:,:)./Bx_vec;
[I, N, T] = size(Bxi_vec);
pvec = reshape(pvec_tmp(:,subset)',[I,1,T]);
run("../subroutines/sub_engel_OLS.m")
[U_vec,WP_total_vec, WP_FS_vec] = CalMoneyMetric_Unobserved(I_vec, Bx_vec, Bxi_vec, pvec, sigmaM, 0);

MM = reshape(U_vec,[N,T]);
Wfull = reshape(I_vec(:,:,:),[N,T]);
[~, indmax] = max(MM);
[~, indmin] = min(MM);
W = Wfull;

for tau = 1:T
Wmax(1,tau) = Wfull(indmax(tau),tau);
Wmin(1,tau) = Wfull(indmin(tau),tau);
wind = (indmin(tau):indmax(tau));
W(:,tau) = NaN;
W(wind,tau) = Wfull(wind,tau);
end


run("../subroutines/MatlabDataConstruction_Raw.m")

T = size(B_vec,3);
N = size(B_vec,2);
J =1000;
U = NaN(1,N,T);
for t = 1:T
    F_UI = interpNaN(log(W(:,t)),log(MM(:,t)));
    U(:,:,t) = exp(F_UI(log(I_vec(:,:,t))));
end

%U_vec = reshape(U,[1,N*T]);
prctiles = prctile(reshape(U,[1,N*T]), linspace(0,100,J)); 


belongs_to = zeros(size(U)); 

N = size(U, 2);
T = size(U, 3);
J = size(prctiles, 2);

for n = 1:N
    for t = 1:T
       
        index = find(U(1, n, t) <= prctiles, 1);
        if isempty(index) 
            index = J; 
        end
        belongs_to(1, n, t) = index;
    end
end

A = NaN(1, J, T); 

By1 = B_vec(10,:,:);
By2 = B_vec(17,:,:);

Sy1 = I_vec.*By1;
Sy2 = I_vec.*By2;


Sy1M = NaN(1, J, T); 
Sy2M = NaN(1, J, T);
WM = NaN(1, J, T);
WM2 = NaN(1, J, T); 



for j = 1:J
    for t = 1:T
        idx = belongs_to(1,:,t) == j; 
        if any(idx)          
            Sy1M(1, j, t) = sum(Sy1(1, idx, t));
            Sy2M(1, j, t) = sum(Sy2(1, idx, t)); 
            WM(1, j, t) = sum(I_vec(1, idx, t))./sum(idx);
            WM2(1, j, t) = sum(I_vec(1, idx, t)); 
        end
    end
end

By1M = Sy1M./WM2;
By2M = Sy2M./WM2;


By1M= reshape(By1M, [], 1);
By2M = reshape(By2M, [], 1);
WM = reshape(WM, [], 1);

year_start = 1974;
T = 44;
J = size(By1M, 1) / T;
year = repmat(year_start:year_start+T-1, [J, 1]); %
year = reshape(year, [], 1); %

bin = repmat((1:J)', [T, 1]);
tbl = table(year,bin, By1M, By2M, WM);

writetable(tbl, '../data/processed/mmquasi.csv');





function F_X = interpNaN(X,Y)
        indexNaN = isnan(X);
        X(indexNaN) = [];% Eliminate NaN
        Y(indexNaN) = [];% Eliminate NaN
        [X, indexUN] = unique(X); 
        F_X = griddedInterpolant(X,Y(indexUN),'linear','none');
end
